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A  Formulation  of  Quantum  Statistical  Mechanics  Based  on  the 
Feynman  Path  Centroid  Density.  I.  Equilibrium  Properties 

Jianshu  Cao  and  Gregory  A.  Voth 

Department  of  Chemistry,  University  of  Pennsylvania.  Philadelphia.  Pennsylvania  19104- 
6323 

ABSTRACT 

A  formulation  of  quantum  statistical  mechanics  is  presented  in  which  the  Feynman 
path  centroid  density  in  Feynman  path  integration  is  recast  as  the  central  statistical  dis¬ 
tribution  used  to  average  equilibrium  and  dynamical  quantities.  In  this  formulation,  the 
path  integral  centroid  density  occupies  the  same  role  as  the  Boltzmarm  density  in  classical 
statistical  mechanics.  Therefore,  the  statistical  ensemble  of  imaginary  time  path  centroid 
configurations  provides  the  distribution  which  is  used  to  average  the  appropriately  formu¬ 
lated  effective  operators  and  imaginary  time  correlation  functions.  An  accurate  renormal¬ 
ized  diagrammatic  perturbation  theory  for  the  centroid  density  and  centroid-constrained 
imaginary  time  propagator  will  also  be  described  with  a  particular  emphasis  given  to  the 
mathematical  advantages  arising  from  the  centroid-based  formulation.  The  present  paper 
is  concerned  with  the  calculation  of  equilibrimn  properties  from  the  centroid  perspective, 
while  the  companion  paper  describes  a  centroid-based  formalism  for  calculating  dynamical 
time  correlation  functions. 
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I.  Introduction 

The  Feynman  path  integrad  representation  of  quantum  mechanics.^  and  in  particular 
of  quantum  statistical  mechanics.’-^  provides  a  powerful  method  to  calculate  and  visualize 
natural  processes.  For  example,  path  integrals  are  particularly  useful  in  describing  the 
quantum  mechanics  of  an  equilibrium  system  because  the  canonical  distribution  for  a 
single  quantum  particle  in  the  path  integral  picture  becomes  isomorphic  with  that  for  a 
classical  “ring  polymer"’  of  quasiparticles^’'^  (cf.  Fig.  1).  In  the  discretized  path  integral 
representation,  the  partition  function  for  a  quantum  particle  is  given  by  the  expression^’"* 

^  (2^^  j  J  dqi---  J  dqp  exp[-dVp{(i)]  .  (1.1) 

where  m  is  the  particle  mass,  (3  equals  l/kpT,  and  the  isomorphic  quasiclassical  polymer 
potential  Vp  is  given  by 

Kp(q)  = 

In  this  discrete  representation,  the  coordinates  {^i}  describe  the  positions  of  the  classi¬ 
cal  quasiparticles  and  have  the  cyclic  property  such  that  qi+p  =  qi  (cf.  Fig.  1).  Each 
quasiparticle  is  harmonically  bound  to  its  two  nearest-neighbors,  and  it  “feels"’  the  inter¬ 
action  potential  through  the  term  V{qi)/P.  In  numerical  apphcations,  a  finite  value  of 
the  discretization  parameter  P  is  used  which  is  large  enough  so  that  suitable  numerical 
convergence  is  obtained.^  While  the  above  equations  have  been  written  for  a  single  quan- 
ttun  particle  in  one-dimension,  a  generalization  to  multiple  particles  and/or  dimensions 
is  straightforward.  In  the  discretized  representation,  the  path  integral  formalism  has 
allowed  for  the  numerical  simulation  of  highly  non-trivial  quantum  systems  using  path 
integral  Monte  Carlo  (PIMC)  or  Molecular  Dynamics  techniques.^ 

In  the  fuUy  quantum  (P  — +  00)  limit,  the  path  integral  for,  e.g.,  the  partition  function 
lo  expressed  as  fimctional  integral  over  all  possible  paths  ^(r)  such  that^ 

Z  =  J  ■■■  j  Vir)  exp{-S[q{T)]/h]}  (1.3) 
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where  the  exponential  weighting  of  the  paths  is  determined  by  the  imaginary  time  action 
functional” 


/•hB 

•^[^('r)]  =  dr  ^  — q{Tf  +  V[q{r)]^  .  (1.4) 

Analytically,  the  path  integral  method  has  proven  to  be  particularly  advantageous  in  the 
theoretical  analysis  of  several  condensed  matter  problems  such  as  the  polaron^  and  the 
solvated  electron.® 

One  of  the  many  interesting  ideas  suggested  by  Feynman  in  his  formulation  and  ap¬ 
plication  of  path  integrals  is  the  notion  of  the  path  centroid  variable,^  denoted  here  by 
the  symbol  qo-  The  centroid  is  the  imagintury  time  average  of  a  particular  closed  Feynman 
path  q{T)  which,  in  turn,  is  simply  the  zero- frequency  Fourier  mode  of  that  path,  i.e.. 


In  the  discretized  path  integral  picture  (i.e.,  for  finite  P),  the  path  centroid  variable  is 
equivalent  to  the  center-of-mass  of  the  isomorphic  polymer  of  classical  quasiparticles  (cf. 
Fig.  1)  such  that 


1  P 

9o  =  •  (1-6) 

1=1 

(It  should  be  noted  that  in  some  previous  publications  the  centroid  variable  has  been 
denoted  by  qo  rather  than  ^o-  The  tilde  notation  will  be  reserved  in  the  present  paper  for 
another  use.) 

Feynman  noted  that  a  quantum  mechanical  “centroid  density”  Pc{qc)  can  be  defined 
for  the  path  centroid  variable  which  is  the  path  sum  over  all  paths  with  their  centroids 
located  at  some  point  in  space  denoted  by  qc-  Specifically,  the  formal  expression  for  the 
centroid  density  is  given  by 


Pc{qc) 


J  ■■■  J  '^q{r)5{qc  -  qo)  exp  {-S{q{T)]/h} 


(1.7) 
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The  quantum  partition  fimction  in  Eq.  (1.3)  is  then  obtained  by  the  integration  of  the 
centroid  density  over  all  possible  positions  of  the  centroid,  i.e., 

Z  =  y  dqcpciqc)  ■  (1.8) 

Some  caution  in  order  when  one  uses  the  centroid  density  because  it  is  distinctly  different 
from  the  coordinate  (or  particle)  density,  p{q),  given  by 

PiQ)  =  =  J  ■■■  j  Vq{T)6[q- qiO)]  exp{-S[qiT)]/h}  ,  (1.9) 

where  g(0)  is  the  beginning  and  end  point  of  the  cyclic  quantum  path  ^(t)  and  is  definitely 
not  the  centroid  variable  in  Eq.  (1.5).  Thus,  the  particle  density  function  is  the  diagonal 
element  of  equiUbrium  density  matrix  in  the  coordinate  representation,  while  the  centroid 
density  does  not  have  a  similar  physical  interpretation.  Nevertheless,  the  integration  over 
either  density  yields  the  quantum  partition  function. 

Feynman  used  the  definition  of  the  centroid  density  along  with  a  simple  approximation 
for  the  action  fimctional  in  Eq.  (1.4)  to  derive  an  expression  for  a  quasiclassical  partition 
function.^  The  latter  function  is  expressed  as  an  integration  over  an  effective  Boltzmann 
factor  which,  in  turn,  depends  on  a  variational  effective  potential  determined  at  each  value 
of  the  centroid  variable  with  the  help  of  an  appropriate  centroid  density  formulation® 
of  the  Gibbs-Bogoliubov  variational  principle.  In  subsequent  work,  several  authors^" 
have  improved  upon  the  Feynman’s  original  approach  by  using  a  more  physically  accurate 
variational  harmonic  reference  system  to  describe  the  imaginary  time  path  fluctuations 
around  the  centroid  variable.  The  effective  harmonic  frequency  and  effective  potential 
are  again  determined  at  each  value  of  the  centroid  variable,  resulting  in  an  approximate 
expression  for  the  centroid  density  in  Eq.  (1.7).  An  approximate  variational  partition 
function®"  can  then  determined  in  such  a  theory  by  virtue  of  Eq.  (1.8). 

One  conclusion  that  can  be  reached  from  the  aforementioned  work  of  Feynman^  and 
others,®" as  well  as  from  the  apparent  utility  of  the  centroid  density-based  formulation 
of  quantum  transition  state  theory, is  that  the  centroid  density  is  a  particularly  useful 
quantity  about  which  to  develop  approximate,  but  highly  accurate,  quantum  mechanical 
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expressions  and  to  probe  the  quantum-classical  correspondence  principle  in  statisticcd  me¬ 
chanics.  It  is  in  this  spirit  that  a  general  centroid  density-based  formulation  of  quantum 
Boltzmann  statistical  mechanics  is  presented  in  the  present  paper.  This  formulation  is 
based  on  a  single  important  notion:  In  any  such  theory  the  centroid  density  should  occupy 
the  same  role  as  the  Boltzmaxm  density  in  classical  statistical  mechanics.  That  is.  the 
rules  for  calculating  operator  averages  and  imaginary  time  correlation  functions  should  be 
reformulated  so  that  the  final  result  of  such  a  calculation  is  obtained  by  an  integration 
over  the  centroid  density  (i.e.,  by  a  statistical  weighting  with  the  centroid  density).  It  wiU 
be  shown  that  such  a  formulation  is  not  completely  trivial  due  to  the  mathematical  differ¬ 
ences  between  the  centroid  approach  and  the  usual  rules  of  quantum  statistical  mechanics. 
Nevertheless,  the  resulting  expressions  are  quite  interesting  both  in  terms  of  the  classical- 
quantum  correspondence  principle  and  from  the  analytical  point  of  view.  One  result  from 
this  approach  is  that  the  centroid  formulation  seems  to  actually  simplify  the  calculation  of 
some  quantum  mechanical  quantities  since  the  imaginary  time  path  fluctuations  have  been 
“pre-averaged”  in  the  centroid  density  expression.  This  statement  is  particularly  true  in 
the  case  of  quantiun  dynamics  (i.e.,  time  correlation  functions)  which  are  discussed  in  the 
companion  paper^^  and  in  Ref.  14. 

Another  significant  feature  of  the  centroid  density-based  formulation  of  statistical 
mechanics  is  that  by  concentrating  on  the  centroid  density  as  the  central  statistical  dis¬ 
tribution.  a  formal  diagrammatic  theory  for  the  centroid  density  can  be  employed  which 
turns  out  to  be  simplified  from  the  point  of  view  of  the  relevant  diagram  topologies.  The 
diagrammatic  expansion  is  therefore  particularly  amenable  to  powerful  renormalization 
techniques.  The  diagrammatic  theory  also  draws  the  formal  connection  between  various 
variational  expressions  for  the  centroid  density  which  have  been  derived  by  others^" and 
specific  diagram  renormalization  strategies.  A  systematic  approach  to  improve  upon  the 
result  of  the  variational  centroid  density  theory  can  then  be  developed.  A  considerable 
amount  of  time  will  be  devoted  to  the  diagrammatic  methods  in  the  present  article  due  to 
the  central  practical  and  formal  importance  of  the  centroid  density  in  the  theory. 

The  sections  of  this  paper  are  organized  as  follows:  In  Sec.  II,  the  basic  equations 
are  derived  for  operator  averages  and  imaginciry  time  correlation  functions  in  the  centroid 
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density-based  formulation  of  quantum  statistical  mechanics.  The  diagrammatic  methods 
for  the  centroid  density  and  related  quantities  are  then  discussed  in  Sec.  Ill,  while  apph- 
cations  are  presented  in  Sec.  IV.  Concluding  remarks  are  given  in  Sec.  V. 

II.  General  Formalism 

In  this  section,  the  general  formulation  of  quantum  Boltzmann  statistical  mechanics 
in  terms  of  the  centroid  density  is  discussed.  The  term  "general  formulation’  means  that 
the  rules  are  presented  here  with  which  to  calculate  operator  averages  and  imaginary  time 
correlation  functions  in  the  centroid  density  perspective.  The  emphasis  will  not  be  on 
new  numerical  procedures  which  might  or  might  not  have  advantages  over  the  already  well 
established  and  successful  numerical  path  integral  techniques  (see,  e.g.,  Ref.  3).  Rather, 
the  theoretical  development  is  more  of  a  conceptual  one  which  is  intended  to  explore  the 
quantum-classical  correspondence  in  statistical  mechanics.  Moreover,  the  rules  outhned 
below  for  calculating  averages  and  correlation  functions  are  formulated  on  the  assumption 
that  a  good  analytic  expression  for  the  centroid  density  can  be  obtained  (cf.  Sec.  Ill), 
thereby  bypassing  the  need  for  a  fully  numerical  approach.  Future  reseeirch  will  be  devoted 
to  an  exploration  of  any  numerical  advantages  inherent  in  the  centroid-based  imaginary 
time  path  integral  formulation.  In  the  companion  paper. it  will  indeed  be  shown  that 
the  centroid-based  approach  provides  a  quite  promising  computational  approach  for  the 
determination  of  quantum  real  time  correlation  functions.^'* 

A.  Imaginary  Time  Position  Correlation  Functions 

At  this  point,  it  is  advantageous  to  introduce  the  notion  of  a  centroid-constrained 
imaginary  time  propagator,  i.e.,  the  correlation  fimction  of  quantum  path  fluctuations 
with  respect  to  the  position  of  the  centroid  variable  9o  =  <7c-  This  correlation  function  is 
defined  as 


X  /•••  /  T>q{r)8{qc-qQ){q{T)-qQ){q{Q)- go)  exv{-S[q{r)]/n) 

J  ■■■  j  Vq(T)6(qc-qQ)exp{-S[q{T)]/n} 

As  the  centroids  of  the  paths  q'(r)  in  this  correlation  function  axe  constrained  to  be  at  qc, 
the  paths  can  be  rewritten  as  q{T)  =  qc  +  Q{T),  where  g(r)  is  the  quantum  path  fluctuation 
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with  respect  to  the  centroid.  A  Fourier  decomposition  of  the  paths  q(r)  can  now  be 
introduced  such  that 


g(r)  =  (2.2) 

n^O 

where  the  smnmation  is  over  all  positive  and  negative  integers  except  for  n  =  0,  and  Qn  is 
Euclidean  frequency  defined  by  fin  =  2-Kn/h0. 

The  correlation  function  in  Eq.  (2.1)  differs  from  the  usual  Euclidean  position  cor¬ 
relation  function  C(r)  because  only  the  paths  with  centroids  at  qc  contribute  to  the 
centroid-constrained  propagator  Cc(t,(Jc).  However,  one  can  obtain  C(t)  by  averaging 
the  centroid-constrained  propagator  over  the  normalized  centroid  density  pc{qc)/Z,  i.e.. 


C(t)  =  (4(t)5(0))  =  (C,(r, ,  (2.3) 

or,  equivalently, 

l(Wr)-9(0;i=)  -  C(0)-C(r)  =  {Cc(0, (Cc(r,gJ)„.  .  (2.4) 

It  should  be  noted  that  the  imaginary  time  function  in  Eq.  (2.4)  provides  a  measure  of 
the  localization  of  quantum  particles  in  condensed  media.^’'^'®  From  this  point  onward, 
the  notation  (•  •  •)p<,  denotes  an  averaging  by  integrating  some  centroid-dependent  function 
over  the  centroid  position  qc  weighted  by  the  normaUzed  centroid  density  Pc{(lc)IZ. 

A  general  method  to  obtain  the  correlation  fimction  Cc{T,qc)  is  through  functional 
differentiation  of  a  generating  functional.  In  order  to  formulate  such  a  fimctional,  an 
imaginary  time-dependent  force  /(r)  is  introduced  into  the  action  functional  of  the  centroid 
density  [Eq.  (1.7)]  such  that  the  potential  V^[g(r)]  is  replaced  by  V^[9('r)]  -t-  /(T)g(r).  The 
centroid-constrained  propagator  is  then  given  by 


C,(r,,,)  =  Urn  1 

/-o  Pc  6f{T)6f{0) 


(2.5) 


Here,  the  centroid  density  pc  is  imderstood  to  be  a  functional  of  the  extra  time-dependent 
force  /(r)  and  is  therefore  the  generating  functionaT^  for  the  centroid-constrained  corre¬ 


lation  function. 
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B.  Operator  Averages 

In  the  normal  path  integral  perspective.--^  the  canonical  average  of  an  operator  -4  is 
given  by  the  expression 


(^4)  =  Z  ^  J  ■  ■  ■  J 'Dq{T)A[q{Q)]  exp{-S[q{r)]/h}  .  (2.6) 

Due  to  the  cyclic  invariance  of  the  trace,  the  operator  can  be  evaluated  at  any  point  along 
the  cyclic  imaginary  time  path  g(r).  This  situation  is  depicted  in  Fig.  2  where  the  centroid 
variable,  which  is  off  the  ring,  is  also  shown.  The  challenge  therefore  is  to  reformulate  the 
rules  for  calculating  the  operator  average  so  that  the  centroid  density  can  be  used  in  the 
statistical  averaging. In  order  to  do  this,  the  average  in  Eq.  (2.6)  is  first  re-expressed  as 


/•••/  Vq{r)6{q^- qo)A[q{r)]e 
f  dqcPciQc) 


The  operator  A  is  then  represented  in  Fourier  space,  i.e.. 


(2.7) 


A{<l)  =  J  ,  (2.8) 

where  A{k)  is  the  Fourier  transform  of  the  operator.  Equation  (2.6)  can  then  be  written 
as 


{A)  =  Z  ^  J  dqcPciqc)  /  ^  ' 


,ikqc 


/■■■/  VqiT)  e^^^^OQ-S[qc+q(r)]/n 


r... 


f  r>g(r) 


,  (2.9) 


where  the  notation  for  the  action  function  “5[gc  -h  <?('r)]”  is  understood  to  denote  path 
fluctuations  around  a  fixed  centroid  at  qc-  At  this  point,  a  cumulant  average  of  the  term 
exp[ikq(T)]  in  the  brackets  in  Eq.  (2.9)  can  be  performed  over  the  path  fluctuation  variable 
q(r).  This  average,  for  the  present  purposes,  is  truncated  at  second-order.  (The  variable 
q(T)  is  also  assumed  to  have  symmetric  fluctuations  about  the  centroid.)  After  performing 
the  inverse  Fomrier  transformation  by  integrating  over  k  in  Eq.  (2.7),  one  arrives  at  the 
flnal  result  of  the  analysis  which  is  given  by 
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(A)  =  {AAqr^),, 


(2.10) 


where  the  effective  centroid-dependent  quasiclassical  function  ,4,.  is  given  by 

ArXqc)  =  {A{qc  +  9))cc(0,<jc) 

^  /n  J  /o  T  /  dqA(qc~q)  expl-q-/2Cci0.qc)]  ■  (2.11) 

y27rCc(0.qc)  J 

Here,  the  variable  q  is  obviously  a  Gaussian  '/axiable  with  a  width  factor  C(;(0.gc)-  The 
effective  classical  function  Adqc)  depends  on  the  centroid  variable  qc  and.  in  order  to 
calculate  the  canonical  average  (.4),  is  averaged  in  a  classical-hke  fashion  in  Eq.  (2.10)  with 
the  normalized  centroid  density  pdqd/Z.  Th'^  above  result  has  been  obtained  previously 
by  one  of  us^^  from  the  somewhat  m*  re  specialized  perspective  of  the  Feynman-Hibbs 
theory.'  If  one  has  an  analytic  or  numerical  representation  of  the  centroid  density  Pc{qc)>.  a 
numerical  quadrature  strategy  can  be  employed  for  evaluating  Eq.  (2.10)  for  cases  in  which 
the  integral  in  Eq.  (2.11)  cannot  be  solved  analytically.  For  multidimensional  systems,  a 
numerical  strategy  would  be  to  use  a  Monte  Carlo  procedure  in  which  Eqs.  (2.10)  and 
(2.11)  are  combined  in  a  single  MC  calculation  based  on  the  importance  sampling  function 

VE(9,gc)  ==  pc{qc)exp[-q  y/2Cc(0.  gc)]  ■ 

It  should  be  noted  that  the  effective  centroid-based  representation  of  the  operator 
-4  in  Eq.  (2.11),  while  useful  in  the  context  of  the  equilibrium  averages,  also  occupies  a 
central  role  in  a  semiclassical  algorithm  to  calculate  quantum  time  correlation  functions  in 
the  centroid  perspective.^^ 

C.  General  Imaginary  Time  Correlation  Functions 

A  general  imaginary  time  correlation  function  for  coordinate-dependent  operators  is 
defined  as^ 

CAB{r)  =  (A[g(T)]S[g(0)]) 

=  J  f  B[q{0)]  exp{-S[q{T)]/h}  .  (2.12) 
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where  the  operators  ^4[g(r)]  and  ^[^(O)]  can  be  evaluated  anywhere  along  the  path  such 
that  0  <  r  <  h3.  This  sduation  is  depicted  in  Fig.  3  and.  again,  the  challenge  is  to 
reformulate  these  rules  so  that  the  centroid  density,  which  is  off  the  ring,  can  be  utilized. 
In  order  to  do  this,  one  can  use  the  fact  that  the  correlation  function  C.^b it'i  can  always 
be  expressed  in  terms  of  centroid-constrained  propagator  Cdr.Qc)  in  Eq.  (2.1)  at  the  level 
of  a  second-order  cumulant  approximation.  To  proceed,  one  first  rewrites  Eq.  (2.12)  as 


C.aslT-) 


fdn 

'/■■■/  Vq{r)  6{q,  -  qo)  .4[q(r)j  B[9(0)]  ' 

1 

f  dqcPciQc) 

(2.13) 


The  operators  A  and  B  are  next  written  in  Fourier  space  as  in  Eq.  (2.8).  Then.  Eq.  (2.12) 
becomes 


CAB(r)  =  Z-'  j  J  ^ 

J  ■■■  f  T>g(T)  ] 

^  I?g{r)  e-^l'7‘=+9dr)l/ft  ’  (2-14) 

The  centroid-constrained  average  in  the  bracket  can  be  performed  as  a  cumulant  expansion 
which,  for  the  present  purpose,  is  truncated  at  second  order.  The  result  for  the  bracketed 
term  is  then  given  by 


{exp[ikiq{T)  +  ik2qi0)])c  ^ 

[klCc{Q,qc)  +2kik2Cc{r,qc)  +  klCc{0,qc)]^  ■  (2.15) 

In  order  to  perform  the  integrals  over  ki  and  kz  in  Eq.  (2.14),  a  coordinate  rotation 
is  first  performed  such  that 


Q+{^)  =  [Qi^)  +  <]{0)]/ 
q-{^)  =  [qi^)~q{0)]/V2 


(2.1b) 


and  new  centroid-constrained  position  correlation  functions  are  defined  as 


Cc{r-qc)  =  Cc{Q.qc)  +  Cc{r,qc) 

C~{r.qc)  =  Cc{Q.qc)  -  Cc{T,qc)  ■  (2.17) 

After  performing  the  integrals  in  Fourier  space,  one  arrives  at  the  final  result  for  the 
imaginary  time  correlation  function  in  the  centroid  density  picture 

CAB{r)  =  {pAB{T,qc))p,  ,  (2.18) 

where  the  centroid-dependent  imaginary  time-correlated  operator  product  pAB(T,qc)  is 
defined  in  the  centroid  density  picture  by  the  double-Gaussian  average 


PAB{r,qc')  =  (^4[(7c  +  2  ^^^{q+  +  q-)]  B[qc  +  2  ^^^{q+ -  q.)]^ .  (2.19) 

Here.  q+  and  <?_  are  Gaussian  variables  with  width  factors  C’^ir^qc)  and  C~{T,qc):  re¬ 
spectively. 

The  double  Gaussian  average  form  for  the  centroid-constrained  operator  product  in 
Eq.  (2.19)  is  somewhat  complicated  to  implement,  particularly  if  the  operators  A  ox  B 
are  not  polynomials  or  exponentials  in  the  variable  q.  It  should  also  be  noted  that  if 
the  cumulant  expansion  is  carried  out  beyond  quadratic  order  in  Eq.  (2.15),  an  even  more 
complicated  analytical  form  of  Eq.  (2.19)  would  be  obtained.  However,  it  is  also  important 
to  remember  that  in  an  analytical  calculation  of  the  imaginary  time  correlation  function 
Cab{'t)  in  the  conventional  path  integral  formalism,  both  the  diagonal  and  off-diagonal 
elements  of  the  thermal  density  matrix  are  required.  The  formahsm  embodied  in  Eq.  (2.18) 
requires  only  an  accmrate  value  of  the  centroid  density  pdqc)  3-nd  the  centroid-constrained 
representation  of  the  correlated  operator  product  in  Eq.  (2.19),  albeit  the  latter  expression 
is  somewhat  complicated.  As  in  the  case  of  operator  averages  (Sec.  IIB),  a  numerical 
quadrature  procedure  can  be  employed  if  Eq.  (2.19)  cannot  be  evaluated  analytically.  In 
the  case  of  multidimensional  systems,  a  combined  Monte  Carlo  procedure  for  evaluating 
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Eq.  (2.18)  simultaneously  with  (2.19)  can  be  developed  using  importance  sampling  based 
on  the  centroid  density  and  the  double-Gaussian  distribution  in  Eqs.  (2.18)  and  (2.19), 
respectively. 

III.  Diagreimmatic  Theory 

-4.  Diagrammatic  Representation  of  the  Centroid  Density 

Because  of  its  central  importance  in  the  centroid-based  theory,  a  systematic  study 
is  presented  in  this  section  of  the  perturbation  expansion,  or  equivalently,  the  cumulant 
expansion  of  the  centroid  density.  This  expansion  has  a  one-to-one  correspondence  with 
a  diagrammatic  representation.  It  will  be  shown  that  diagrammtic  classifications  and 
topological  reductions  result  in  the  renormalization  of  the  vertices  and  lines,  and  thus  lead 
to  a  set  of  self-consistent  equations.  Previous  results  based  on  an  optimized  harmonic 
reference  potential^" can  be  readily  derived  and  thus  systematically  improved. 

In  general,  it  will  be  assumed  that  the  Euclidean  action  of  the  reference  system  takes 
a  quadratic  form  in  the  variable  q(T}  such  that 

Sref[q{r)]  = 


where  ctn  defines  the  reference  centroid-constrained  propagator  [cf.  Eq.  (2.1)]  such  that 

a(r)  =  ^a„e-""'  .  (3.2) 

and  Qn  =  0  due  to  the  centroid  constraint.  Two  well-known  quadratic  models  are  the  free- 
particle  reference  system,  where  a“^  =  m/3Q^,  and  the  linear  harmonic  oscillator  (LHO) 
reference  system,  where  =  ml3{Q^  with  u  being  the  intrinsic  LHO  frequency. 

With  a  solvable  reference  system  in  hand,  one  can  express  the  centroid  density  as 

PciQc)  =  Pc,ref{gc)  {exp  (3.3) 
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where  AV^  is  the  imaginary  time  average 

_  1 

^  MJ,,  + 9(^)1  :  (3-4) 

and  AV"  =  V —Vref  is  the  deviation  of  the  real  potential  V  from  the  reference  potential  VVe/- 
The  symbol  (•  ■  ■)c,ref  in  Eq.  (3.3)  denotes  a  centroid-constrained  path  integral  average  in 
the  reference  system.  As  the  centroid-constrained  propagator  of  the  reference  potential 
Q(r)  uniquely  defines  an  infinite  set  of  Gaussiaii  averages  over  the  Fourier  modes  {qn}, 
one  can  equivalently  denote  the  centroid-constrained  average  in  the  reference  system  with 
the  symbol  (•  ■  ■)q. 

The  first  step  in  the  development  of  a  diagrammatic  theory  for  the  centroid  density 
is  to  Taylor  expand  the  average  in  Eq.  (3.3),  i.e., 


_  °°  I  /{  I  ' 

{exp{-0EV))a  =  dTAV{qc  +  qir))\ 


where  AV{k,qc)  is  the  spatial  Fourier  transformation  of  diflFerence  potential  AV(gc  +  q) 
with  respect  to  the  variable  q.  It  should  be  noted  that  the  subscript  n  in  Eq.  (3.5)  differs 
from  that  used  in  the  Fourier  expansion  in  Eq.  (3.2).  Since  q{Ti)  in  the  reference  system 
is  a  Gaussian  variable  with  zero  mean  at  any  imaginary  time  Xj,  the  cumiilant  expansion 
of  a  linear  combination  of  those  variables  trvmcates  at  second  order,  giving 


\  S  A:2q(0)  +  ^  ^  k^kja{Ti  -  r,)  i  ,  (3.6) 

i=l  0^3  J  J 

in  which  a(r)  is  defined  by  Eq.  (3.2).  There  is  no  lineeir  term  in  Eq.  (3.6)  because  (g)  =  0 
according  to  the  definition  of  the  variable  g(r)  being  the  path  fluctuation  with  respect  to 
the  centroid.  This  property  greatly  simplifies  the  cumulant  expansion  and  diagrammatic 
analysis.  By  substituting  the  Taylor  expansion  of  Eq.  (3.6)  into  Eq.  (3.5),  and  transforming 
back  into  coordinate  space,  one  arrives  at  the  result 


U=1 


,ikiq{Ti) 


=  exp  < 
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(exp  (— /3AV))^gj  =  EE 

n=0  m=0 


(-ir 

n\m\ 


"  />ft/3  j_,  1  ^  1 

n  /  ~r  ~ 


AV  ,  {3.7a) 


where 


AV  s  nAV'(9,  +  ?(Ti)l|^.^„  (3.76) 

1=1 

and  the  partial  derivative  symbol  is  defined  to  be  di  =  d/dq{Ti)  appfied  at  timesUce  t^. 
The  imaginary  time  integrals  in  Eq.  (3.7a)  are  understood  to  be  integrations  over  any  and 
all  imaginary  time  slices  which  appear  in  the  expansion.  After  the  derivatives  are  taken, 
the  potential  difference  terms  are  evaluated  at  the  position  of  the  centroid  Qc- 

In  principle,  the  above  perturbation  seri^  gives  a  complete  description  of  the  quantum 
mechanical  centoid  density,  and  the  expansion  can  be  ceirried  out  to  any  order  according  to 
the  accuracy  required  for  a  specific  problem.  However,  it  is  conceivable  that  in  higher  orders 
the  expansion  terms  become  increasingly  complicated  and  a  low-order  calculation  will  not 
provide  a  reasonable  physical  picture  of  a  quantum  system.  Therefore,  a  diagrammatic 
representation  of  Eq.  (3.7)  is  introduced  here  which  provides  a  powerful  way  to  analyze 
the  pertvubation  expression  and  to  visualize  the  analytical  expressions. 

A  careful  examination  of  the  series  in  Eq.  (3.7)  reveals  the  basic  composition  of  the 
expansion  terms.  That  is,  all  diagrams  consist  of  two  basic  elements,  vertices  emd  lines. 
Each  vertex  is  designated  by  a  Euclidean  time  Tj  to  be  integrated  from  Tj  =  0  to  Tj  =  ^/3, 
and  the  potential,  or  its  derivatives,  axe  evaluated  at  the  position  of  the  centroid.  Each  line 
connecting  two  vertices  at  times  Tj  and  Tj  is  designated  as  a  reference  centroid-constrained 
propagator  Q:(Tt  —Tj).  Whenever  a  Une  connects  to  a  vertex,  a  spatial  derivative  is  applied  to 
the  potential  so  that  the  order  of  the  derivative  is  equal  to  the  number  of  lines  that  connect 
to  the  vertex.  A  negative  sign  is  assigned  to  each  vertex.  The  value  of  a  diagram  is  the 
product  of  all  the  composing  elements  which  multiply  a  symmetry  coefficient  determined 
by  the  topological  structure  of  the  diagram. 
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With  above  definitions  in  hand,  one  can  establish  a  one-to-one  correspondence  between 
each  distinct  perturbation  term  and  each  diagram.  The  expansion  series  of  Eq.  (3.7)  is 
then  the  collection  of  all  topologically  different  diagrams  and  all  possible  combinations. 
The  following  diagrams  are  some  examples: 


o  ^  r\ 

_  1  +  o  +  +  O  +  (  + 

{exp{  —  3AV))a  =  o  O  Xj 


O  +  +  •  ■  •  . 


(3.8) 


A  well-known  graph  theorem^®  states  that  an  infinite  series  of  all  possible  topologi¬ 
cally  different  diagraims  and  their  combinations  is  equal  to  the  exponential  of  all  possible 
topologically  different  connected  diagrams.  For  a  connected  diagram,  any  two  vertices  axe 
linked  to  each  other  by  at  least  one  line  or  one  path  of  fines.  Therefore,  one  can  to  express 
the  centroid  density  as 


PciQc)  =  Pc,ref{qc)  exp{T) 


(3.9) 


where  JE  is  given  by 

T 

The  underlying  diagranunatic  techniques  employed  here  are  not  new,  having  appeared 
many  times  in  the  literatme  such  as  in  the  Meyer  cluster  expansion^^  and  the  Fejoiman 
diagrams.^®  However,  in  order  to  apply  these  techniques  to  a  specific  problem  one  must 
search  for  the  proper  diagrammatic  representation.  In  the  case  of  the  centroid  density,  the 
present  approach  is  not  restricted  by  the  functional  form  of  the  potential  and  proves  to  be 
particularly  advantageous  in  analytical  studies.  It  should  also  be  noted  that  the  topological 
reduction  performed  in  the  present  case  is  equivalent  to  a  diagrammatic  representation  of 
the  cumulant  expansion.  The  cumulant  expansion  itself  becomes  tedious  in  high  orders 
and  there  are  a  large  number  of  cancellations  in  the  present  approach  not  explicilty  given 
in  the  cumulant  relations.  It  therefore  proves  to  be  much  easier  to  keep  track  of  higher 
order  terms  in  the  diagrammic  representation. 


0  ♦  oo  -  o 


(3.10) 
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As  pointed  out  earlier,  all  the  diagrams  of  T  are  closed  due  to  the  fact  that  q  =  0 
(i.e..  the  centroid  constraint).  For  example,  there  exists  no  class  of  diagrams  D  with  single 
lines  hanging  outside  of  the  main  diagram,  such  as 

^  ^  t  +  >>^0=  •  (3,U) 

These  kinds  of  diagrams  contribute,  e.g.,  to  the  expansion  for  the  usual  quantum  den¬ 
sity  {q\exp{—l3H)\q).  This  feature  of  Eq.  (3.10)  simplifies  the  analysis  enormously  and 
makes  the  centroid  density  perspective  a  prefereable  choice  for  the  analytical  calculation 
of  equilibrium  averages  and  imaginary  time  correlation  functions  in  quantum  statistical 
mechanics.  It  remains  to  be  seen,  however,  whether  the  centroid  perspective  can  offer 
any  numerical  advantage  over  standard  numerical  path  integral  approaches  for  calculating 
equihbrium  properties.^  This  issue  will  be  explored  in  future  research. 

B.  Renormalization  of  the  Vertices 

The  diagrammatic  representation  of  the  centroid  density  enhances  our  intuition  and 
capability  in  an  approximate  evaluation  of  the  full  perturbation  series.  Instead  of  a  tedious 
term-by-term  calculation,  one  can  focus  on  a  class  of  diagrams  with  the  same  topological 
characteristics.  The  sum  of  such  a  class  often  results  in  a  compact  analytical  expression 
which  includes  infinite  terms  in  the  summation.  A  very  useful  technique  in  such  cases 
is  the  renormalization  of  diagrams. This  procedure  is  first  appUed  here  to  the  vertices 
to  obtain  an  effective  potential  and  thereby  an  accurate  approximation  to  the  centroid 
density. 

The  first  set  of  diagrams  to  be  studied  contain  only  one  vertex,  i.e., 

o  =  o  +  Q  +  OO  +  ••• 

=  -/3AF-i/3AK(2)^(0)_^^^^(4)^2(o)  +  ...  ^  (3.12) 

where  the  superscripts  “(z)”  denote  the  order  of  the  spatial  derivative  anJ  all  terms  AV 
and  AV^*^  are  evaluated  at  the  centroid  :  '’sition  qc-  The  different  terms  in  Eq.  (3.12) 
correspond  to  the  corrections  due  to  local  quantum  path  fluctuations  q{T).  Summing  up 
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the  series,  one  obtains  a  closed  form  expression  for  Eq.  (3.12)  in  the  form  of  a  Gaussian 
average  over  a  single  variable  with  a  Gaussian  width  factor  q(0)  such  that 


^  ~1^==7^  f  +  exp[-q-/2a{0)]  (3.13) 

\/ 27ra(U)  J 

This  form  of  the  effective  potential  incorporates  to  a  certain  degree  quantum  effects  and 
the  anhaxmonicity  of  the  potential.  In  the  case  of  the  free  particle  reference  system,  the 
centroid-constrained  propagator  is  given  by 

atp{r)  =  ^13(1  -  2uf  -  1)  (3.14) 

where  u  =  t/TiQ  and  =  h}0/\2m.  Substitution  of  this  expression  for  q;(0)  in  Eq.  (3.13), 
and  using  fact  that  for  the  free  particle  reference  system  Vref{q)  =  0,  one  recovers  the 
well-known  Feynman-Hibbs  quasiclassical  theory^  for  centroid-dependent  effective  classi¬ 
cal  potential  from  Eq.  (3.13).  Of  course,  any  general  quadratic  reference  system  for  the 
propagator  a{r)  can  be  used  in  the  present  theory  and  will  lead  to  greater  accuracy  (see 
below). 

The  next  order  diagram  to  be  considered  is  a  ring  diagram  with  two  vertices  and  two 
lines,  i.e., 

dn  /  dT2a^{Ti-T2)  .  (3.15) 

Ah‘‘  Jo  Jo 

As  an  example,  substituting  Eq.  (3.14)  into  Eq.  (3.15)  yields  the  first  correction  term  to 
the  Feynman-Hibbs  approximation 

1  rhff  rh0  /?2\4 

_(^y(2))2^  dn  dT2a}p{n-T2)  =  .  (3.16) 

which  is  exactly  what  one  would  expect  from  the  cumuleint  expansion  [see,  e.g.,  Eq.  (2.28) 
in  Ref.  10].  Again,  any  quadratic  reference  propagator  q:(t)  could  be  used  here  instead  of 
the  free  particle  one. 

A  further  correction  is  to  add  local  quantum  fluctuations  to  the  diagrams  with  two 
vertices,  the  latter  diagrams  being  given  by 
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<0  =  <0 -00*000 


(3.17) 


This  procedure  is  equivalent  to  replacing  the  potential  AV  in  Eq.  (3.15)  with  the  effective 
potential  {AV)a  of  Eq.  (3.13).  This  correction  can  be  introduced  in  the  same  fashion  for 
all  higher-order  diagrams.  However,  before  doing  this  it  is  advantageous  to  include  more 
diagreims  in  the  vertex  corrections.  The  set  of  local  fluctuation  diagrams,  Eq.  (3.12),  is  only 
the  simplest  correction.  One  can  extend  the  analysis  to  incorporate  all  ring  corrections, 
given  by 


1  ^fi/3  fh0 

-f-  /  dTi  /  dT2a{Ti-T2)a{r2-Ti){AV^'^'>)a 

2n  jQ  Jo 

+  ••• 


(3.18) 


By  expressing  the  convolution  integrals  in  the  powers  of  ocn  >  one  can  obtain  a  simple  closed 
form  for  the  ring  diagrams,  given  by 


-0 


E, 


a: 


n#0 


21+0(AV(^% 


AV  =  -0VAV  , 


(3.19) 


where  V  should  be  understood  as  an  operator.  Furthermore,  one  can  include  diagrams 
with  multiple  rings  hanging  on  the  same  vertex  which  symbolically  leads  to 


o  =  o  +  (^  +  +  ■■■ 

=  -0(l+V-i-^V  +  ^V^  +  ---)AV  .  (3.20) 

Equation  (3.20)  now  deflnes  an  effective  potential  which  has  the  same  form  as  Eq.  (3.13), 
except  the  Gaussian  width  factor  is  now  given  by 
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(3.21) 


_ _ 

This  is  a  particularly  useful  form  of  the  equation  for  the  effective  potential  because  it 
involves  only  the  second  derivative  of  the  potential  averaged  about  the  centroid  using  the 
unoptimized  reference  propagator  width  factor  a. 

As  has  already  been  suggested.  AV  in  Eq.  (3.21)  can  be  replaced  by  an  effective  po¬ 
tential  difference  AV.  The  procedure  here  is  called  the  renormaUzation.  Analyticedly,  the 
renormalized  element  appears  as  an  unknown  variable  to  be  solved  from  self-consistent 
equations  which  relate  the  given  unrenormalized  quantities  to  the  renormalized  ones.  Usu¬ 
ally,  it  is  straightforward  to  discover  such  relationship  by  observing  the  topologj-  of  dia¬ 
grams.  As  a  result  of  the  renormaUzation  of  the  vertices,  one  arrives  at  the  expression 
for  the  renormaUzed  potential  difference  AV  by  replacing  the  Gaussian  width  factor  a  in 
the  averaging  in  the  right-hand-side  of  Eq.  (3.21)  by  the  effective  width  factor  a.  Note 
that  this  equation  must  now  be  solved  self-consistently.  The  resulting  expression  for  AV 
is  given  by 


AU  =  {AV{qc  +  q))a  ,  (3-22) 

with 


-  I3AV  =  • 

and  (•  •  •)q  denotes  a  Gaussian  average  with  a  width  factor 


(3.23) 


a 


=  E 

n^O 


OCn 

1  -b  /?AU(2)a„ 


(3.24) 


The  notation  means  that  the  second  derivative  of  AV{qc+q)  is  taken  with  respect 

to  q.  RenormaUzed  quantities  wiU  be  denoted  hereafter  with  an  overbar.  Diagrammati- 
cally,  a  black  vertex  stands  for  a  renormaUzed  potential  [cf.  Eq.  (3.23)]  while  a  bold  line 
stands  for  a  renormaUzed  centroid-constrained  propagator.  Again,  the  underlying  logic  of 
renormaUzation  is  not  new  and  has  been  used  many  times  in.  e.g..  Green’s  function  theory. 


-  19  - 


mean  field  theory,  etc.  To  our  knowledge,  it  has  not  been  applied  in  a  general  way  to  treat 
the  centroid  density.  A  multidimensional  generalization  is  discussed  in  the  Appendix. 

The  seemingly  complicated  equations  given  above  actually  have  a  relatively  simple 
interpretation.  By  substituting  the  centroid-constrained  propagator  for  the  LHO  into 
Eq.  (3.24),  and  using  a  general  LHO  frequency  a;  such  that  AV^  =  V  —  \mu:-q^  and 
=  m3{Q‘^-r'jj-),  the  renormalized  LHO  frequency  ui  can  be  specified  firom  the  definition 


=  {V^^Hqc  +  q))a  ■  (3.25) 

The  value  of  the  renormalized  LHO  frequency  is  determined  from  the  solution  to  the  above 
transcendental  equation,  where  the  renormalized  Gaussian  width  factor  in  Eqs.  (3.24)  and 
(3.25)  is  now  given  by 


(3.26) 


It  should  be  noted  that  these  equations  are  to  be  solved  for  each  position  of  the  cen¬ 
troid  qc.  The  frequency  in  Eq.  (3.25)  is  exactly  the  effective  frequency  obtained  for  the 
optimized  LHO  reference  system  using  the  path  integral  centroid  density  version  of  the 
Gibbs- Bogoliubov  variational  method.®  Originally  suggested  by  Feynman.®'®  this  varia¬ 
tional  theory  has  also  been  employed  in  the  evaluation  of  quantum  rate  constants  and 
to  improve  the  convergence  of  path  integral  Monte  Carlo  simulations.^®  The  present  deriva¬ 
tion,  however,  does  not  depend  on  a  specific  reference  system  as  long  as  it  is  quadratic.  In 
addition,  Eqs.  (3.25)  and  (3.26)  are  not  derived  from  any  variational  principle,  but  are  in¬ 
stead  the  result  of  diagram  renormalization.  More  importantly,  the  diagrammatic  analysis 
provides  a  way  to  systematically  improve  upon  the  variational  theory. 

In  order  to  improve  upon  the  optimized  LHO  theory  for  the  effective  centroid  potential, 
one  need  to  consider  the  contribution  from  higher  order  diagrams.  One  way  to  accomplish 
this  is  to  diagrammatically  impose  the  condition  (AV’^^^)q  =  0  [cf.  Eqs.  (3.25)  and  (3.26)]. 
This  condition  specifies  that  all  vertices  linked  to  two  lines  will  'vanish,  i.e.. 
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Consequently,  all  diagrams  containing  this  element  will  vanish,  giving  the  result  a  =  a. 
The  leading  corrections  in  the  centroid  density  expansion  in  Eq.  (3.10)  are  then  given  by 


(3.28) 


and 


where  the  centroid-constrained  propagator  for  the  LHO  is  given  by 


d(r)  = 


m0u>'^ 


b/2 


cosh(l  —  2u)b/2  —  1 


(3.29) 


(3.30) 


_sinh(6/2) 

in  which  b  =  hf3uj  and  u  =  r/^/3.  These  two  terms  provide  an  improvement  on  the 
optimized  harmonic  reference  centroid  density  approximation,  giving 


Pc  -  Po,re/exp|  /3[(AC)a  213' (mil-) 

_  ( AS]  (A1/W)2  +  •  •  ■]  1 .  (3.31) 

2!4!  hu  \mujj  ^  Jj  ^  ^ 

Here,  /n(fc)  is  a  dimensionless  coefficient,  defined  by 

which  becomes  a  constant  in  the  limit  of  large  6.  More  corrections  can  be  included  by 
adding  more  diagrams.  This  procedure  will  be  discussed  in  the  next  subsection  in  the 
context  of  renormalization  of  the  centroid-constrained  propagator  (i.e.,  the  lines). 

C.  Renormalization  of  the  Lines 

The  other  essential  element  of  the  diagrams  is  the  line  which  represents  the  centroid- 
constrained  propagator.  The  renormalization  of  vertices  leads  to  the  evaluation  of  centroid 
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density  while  the  renormalization  of  lines  leads  to  the  evaluation  of  Euclidean  centroid- 
constrained  propagator  defined  in  Eq.  (2.1).  This  propagator  can  be  obtained  formally 
from  Eq.  (2.5)  combined  with  Eq.  (3.9). 

Following  the  diagrammatic  analysis  of  Sec.  2.  all  the  leachng  diagrams  for  Cdr.qc) 
are  obtained  from  Eq.  (3.10),  i.e.. 


Qc)  — 


Obviously,  the  above  collection  of  diagrams  represents  all  possible  contributions  to 
the  centroid-constrained  correlation  function.  In  fact,  all  the  decorations  attached  to  the 
intermediate  vertices  can  be  removed  if  the  vertex  is  renormalized.  This  operation  can  be 
achieved  by  replacing  all  the  AK’s  by  AF’s,  giving 

AV-  =  (AVfe+fl)c,(o.,.)  .  (3.34) 

where  the  Gaussian  width  factor  is  now  Cc{0,qc)  instead  of  q(0).  In  the  case  of  fully  renor¬ 
malized  vertices,  Cc(0,  qc)  is  equivalent  to  the  renormalized  reference  centroid-constrained 
propagator  q  of  Eq.  (3.24).  These  two  notations  will  not  be  distinguished  hereafter. 

The  simplest  set  of  lines  is  the  chain  collection,  given  by 

-  =  +  •  +  •  •  +  •  •  •  .  (335) 

where  the  bold  line  stands  for  a  renormalized  line.  This  diagram  can  also  be  expressed  in 
the  compact  form 


which  leads  to  the  following  self-consistent  equation, 

dn  =  an-  PAV^'^^Otnan 


(3.36) 


(3.37) 
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By  noting  that  q  =  (3-34)  and  Eq.  (3.37)  are  the  same  as  the 

optimized  LHO  reference  equations  in  Eqs.  (3.22)  and  (3.24)  derived  in  the  last  section.  In 
fact,  differentiating  a  ring  collection  of  the  generating  functional  always  produces  a  chain 
collection  in  the  correlation  diagrams.^®  The  multidimensional  generalization  of  the  above 
equation  is  described  in  the  Appendix. 

The  next  stage  in  the  analysis  is  to  include  all  the  two-liut^loop  corrections  in  the 
renormalization,  given  by 


=  .(3.38) 

Because  of  the  imaginary  time  convolutions  in  the  expression,  the  analytical  expressions 
for  the  above  diagrams  are  written  in  Fourier  space,  where  a*n  is  the  contribution  from 
the  two-line-loop,  given  by 


Q^n  —  'y  ^  ^n~m^m  ■  (3.39) 

m^O 

Since  Q“  is  a  convolution  expression,  the  self-consistent  equation  for  a  is  not  local  in  Fourier 
space,  and  therefore  we  can  no  longer  seek  a  single  effective  frequency  solution  as  in  Eq. 
(2.15).  As  a  matter  of  fact,  this  analysis  shows  that  the  optimized  LHO  reference  system 
reaches  the  maximum  capacity  for  a  quadratic  potential  to  approximate  an  anharmonic 
potential  and  any  further  corrections  are  beyond  an  effective  frequency  description. 

The  infinite  smnmation  in  Eq.  (4.10)  can  be  carried  out  to  yield  a  closed  equation, 
given  by 


Qfn 


1  + 


(3.40) 


which  can  be  solved  numerically.  It  is  important  to  incorporate  infinite  terms  corresponding 
to  the  same  class  of  diagrams  so  that  at  low  temperature  and  high  aharmonicity  the  self- 
consistent  equation  will  not  diverge. 
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Further  corrections  to  the  propagator  diagrams  will  consist  of  multi- line- loops  and 
their  combinations.  However,  the  expressions  become  increasingly  complicated  and  we 
have  not  been  able  to  reach  a  general  result  for  the  infinite  summation  of  all  the  multi¬ 
line-loops.  Apart  from  the  simple  multi-line-loops,  there  also  exist  a  large  number  of 
irreducible  diagrams  that  are  not  included  in  the  above  renormahzation  scheme.  It  may 
be  possible  that  information  from  some  of  the  lower-order  diagrams  could  help  in  the 
construction  of  an  accurate  renormalization  equation  by  means  of  a  Fade  apprcximant  or 
a  continued  fraction  scheme. 

As  shown  in  Sec.  II.  the  centroid-constrained  Euclidean  correlation  function  is  of 
central  importance  in  studying  equilibrium  properties  via  the  centroid  density  perspective. 
Even  more  importantly,  its  real-time  counterpart  is  essential  in  describing  the  quantum 
dynamics  of  a  system  (cf.  the  companion  paper^^  and  Ref.  14) .  The  real  time  and  imaginary 
time  correlation  functions  are  of  course  related  by  the  analytical  continuation  r  ^  it.  so 
the  det£Liled  study  of  the  Euclidean  centroid-constrained  correlation  function  presented  in 
this  section  may  eventually  help  us  to  understand  the  real  time  behaviour  of  quantum 
systems. 

IV.  Applications 

In  this  section,  the  results  of  calculations  which  probe  the  accuracy  of  the  four  main 
topics  discussed  in  this  paper  are  presented.  These  topics  are:  (1)  the  centroid-based 
formulation  for  calculating  equilibrium  averages  (Sec.  IIB);  (2)  the  centroid-based  for¬ 
mulation  for  calculating  imaginary  time  correlation  functions  (Sec.  IIC):  (3)  the  analytic 
diagrammatic  approach  for  the  calculation  of  the  centroid  density  (Sec.  IIIB);  and  (4)  the 
analytic  diagrammatic  approach  for  calculating  the  imaginary  time  propagator  (Secs.  IIA 
and  IIIC). 

For  all  of  the  above,  the  numerical  calculations  are  based  on  a  completely  non¬ 
quadratic  potential,  given  by 


V(q)  =  ,3  +  ,V2 


(4.1) 


where  the  mass  m  and  h  are  taken  to  be  unity.  All  of  the  analytical  results  are  also 
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ri^'inpared  with  PIMC  simulation  results.  To  achieve  good  convergence,  the  path  integral 
simulations  emoloyed  P  =  100  discretizations  (or  "polymer  "  quasiparticles)  and  MC  10” 
passes.  The  nuntber  of  beads  moved  on  each  trial  was  adjusted  to  yield  an  acceptance  rate 
(d"  oOVy. 

In  Table  I.  the  values  tor  the  equilibrium  average  {q‘)  are  tabulated  for  the  poten¬ 
tial  in  Eq.  (4.1)  as  calculated  by  PIMC  and  by  the  centroid  formulation  in  Eq.  (2.10) 
with  the  optimized  LHO  approximation  to  the  centroid  density  [cf.  Eqs.  (3.22)-(3.26)]. 
Also  tabulated  are  the  analytic  results  obtained  using  the  optimized  LHO  theory  with 
the  higher-order  corrections  in  Eqs.  (3.28)  and  (3.29).  The  results  are  shown  for  various 
temperatures.  Clearly,  the  centroid  approach  is  an  accurate  one.  It  should  also  be  noted 
that  these  calculations  represent  a  stringent  test  of  the  centroid  formulation  because  the 
potential  in  Eq.  (4.1)  contains  no  intrinsic  quadratic  term. 

In  Fig.  4.  the  imaginary  time  correlation  function  {q~ {T)q~ [0])  is  plotted  for  a  temper¬ 
ature  of  =  5  and  as  a  function  of  the  dimensionless  variable  u  =  rjhl3.  The  solid  circles 
depict  the  exact  PIMC  result  while  the  solid  line  is  for  the  optimized  LHO  result  along  with 
the  centroid-based  formulation  of  the  correlation  function  in  Eq.  (2.18).  In  Fig.  5.  simi¬ 
lar  results  are  shown  for  the  correlation  function  (q^(T)q^(O)}.  Again,  the  centroid-based 
formalism  is  in  very  good  agreement  with  the  numerically  exact  result. 

Although  the  calculations  described  above  provide  an  indirect  test  of  the  analytic 
expressions  derived  for  the  centroid  density,  it  is  desireable  to  provide  a  direct  test.  There¬ 
fore,  the  quantum  correction  factor  for  a  Eckart  barrier  has  been  calculated  and  tabulated 
in  Table  H  within  the  context  of  path  integral  quantum  transition  state  theory.  In  the 
context  of  that  theory,  the  quantum  correction  factor  T  to  the  classical  rate  constant  is 
given  by 


r  =  Pc(9*)/Pd(9*)  ,  (4.2) 

where  Pc  and  pd  are  the  centroid  density  and  the  classical  density  at  the  barrier  top  {q  = 
q' ) ,  respectively.  Since  the  centroid  density  is  larger  than  the  classical  density  at  the  barrier 
top,  the  quantum  rate  is  enhanced  by  the  correction  factor  T.  In  their  paper.  Voth  et  al.^^“ 
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evaluated  this  factor  for  an  Eckart  barrier  at  different  temperatures  and  found  the  difference 
between  the  numerical  path  integral  results  and  the  optimized  LHO  approximation  for 
the  centroid  density  increases  as  the  temperature  is  lowered.  Therefore,  the  higher-order 
analytic  corrections  in  Eq.  (3.31)  become  important  in  the  extreme  quantum  limit.  As  the 
Eckart  barrier  is  an  even  potential,  the  cubic  term  in  Eq.  (3.31)  vanishes  at  the  barrier 
top  and  the  leading  correction  is  quartic.  In  Table  II,  the  results  are  hsted  from  the 
optimized  LHO  reference  approximation,  from  the  optimized  LHO  approximation  including 
the  higher  order  corrections,  and  from  the  PIMC  results.  For  the  last  two  entries,  the 
quartic  correction  seems  to  slightly  overestimate  the  PIMC  result,  although  in  either  case 
the  analytic  theory  is  essentially  exact  to  within  the  MC  error  (~  5%)  of  the  numerical 
results. 

A  final  set  of  calculations  was  designed  to  test  the  accuracy  of  the  analytic  expressions 
for  the  centroid-constrained  correlation  function  Cdr^qc)  [Eq-  (2.1)],  as  well  well  as  the 
imaginary  time  position  correlation  function  C{t)  obtained  by  averaging  the  centroid- 
constrained  propagator  over  the  coordinate  space  weighted  by  the  centroid  density.  To 
demonstrate  the  effect  of  the  higher-order  correction  terms  in  Eq.  (3.40),  the  centroid- 
constrained  propagator  was  first  evaluated  using  the  optimized  LHO  approximation  in  Eq. 
(3.37)  as  shown  by  the  dashed  line  in  Fig.  6  for  /?  =  10  and  Qc  =  0.  The  result  obtained 
by  including  the  two-line-loop  correction  in  Eq.  (3.40)  is  shown  by  the  sohd  line,  while  the 
PIMC  results  are  given  by  the  solid  circles.  The  results  are  again  calculated  as  a  function 
of  It  =  r/hp.  Though  the  effective  LHO  approximation  provides  a  good  approximation  to 
the  exact  centroid-constrained  correlation  function,  the  correction  from  the  two-line-loop 
diagram  renormalization  clearly  improves  the  agreement  with  the  numerical  data.  In  fact, 
the  latter  theoretical  prediction  virtually  coincides  with  the  exact  PIMC  result. 

Figure  7  shows  the  imaginary  time  position  correlation  function  C(r)  for  /?  =  5. 
The  effective  frequency  ui  was  solved  self-consistently  from  Eqs.  (3.25)  and  (3.26)  at  each 
centroid  position  Qc  and  then  the  centroid-contrained  optimized  LHO  correlation  function 
was  averaged  with  the  numerical  centroid  density.  As  seen  from  Fig.  7,  good  agreement 
was  obtained  with  the  numerical  result.  Even  better  agreement  was  obtained  from  Eq. 
(3.40)  using  the  two-line-loop  correction  (solid  curve).  Again,  the  deviation  from  the  exact 
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result  in  either  case  is  quite  small. 

V.  Concluding  Remarks 

In  the  present  paper,  the  basic  computational  procedures  of  equilibrium  statistical 
mechanics  have  been  reformulated  in  order  to  cast  the  path  integral  centroid  density  into 
role  of  the  underlying  statistical  distribution  in  quantum  statistical  mechanics.  The  usual 
expressions  for  equihbrium  averages  and  imaginary  time  correlation  functions  have  been 
modified  so  that  the  final  answer  is  obtained  by  averaging  an  effective  centroid-dependent 
function  over  the  centroid  distribution.  These  effective  functions  involve,  in  the  case  of 
an  operator  average,  a  quasiclassical  centroid-dependent  function  or,  in  the  case  of  an 
imaginary  time  correlation  function,  a  quasiclassical  operator  product.  In  each  case,  the 
quasiclassical  centroid-dependent  quantity  is  formulated  as  a  Gaussian  averaged  function 
“broadened”  by  the  intrinsic  quantum  thermal  width  of  the  particle. 

In  addition  to  the  computational  formalism  for  averages  and  correlation  function,  a 
diagrammatic  perturbation  theory  for  the  centroid  density  and  related  quantities  has  been 
formulated  and  analyzed  in  some  detail.  The  connection  between  different  levels  of  dia¬ 
grammatic  summation  and  renormalization  in  the  perturbation  theory  and  the  Feynman- 
Hibbs^  or  variational  theories^"^^  for  the  centroid  density  have  been  identified,  leading 
to  a  systematic  methodology  for  improving  upon  the  latter  two  approaches.  Correspond¬ 
ingly,  specific  correction  factors  have  been  derived  and  excellent  agreement  with  numerical 
calculations  obtained. 

The  primary  motivation  for  the  development  of  the  present  formalism  eirises  from  the 
numerous  appealing  properties  of  the  centroid  density  such  as  its  compelling  analogy"  with 
the  classical  Boltzmaan  density  and  the  topological  simplicity  of  its  imderlying  diagram¬ 
matic  representation.  Another  motivation  is  the  apparently  central  role  occupied  by  the 
centroid  density  in  the  path  integral  quantum  transition  state  theory^^  for  activated  rate 
constants.  These  facts  have  lead  us  to  develop  the  more  general  perspective  presented 
in  the  present  paper.  The  companion  paper^^  explores  the  intriguing  role  played  by  the 
centroid  variable  and  centroid  density  in  dynamical  quantum  time  correlation  functions.^"* 
Applications  of  the  centroid-based  formalism  described  in  the  present  and  companion  pa- 
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pers  will  be  the  subject  of  future  research. 


Acknowledgements 

This  research  was  supported  by  the  National  Science  Foundation  and  the  Office  of 
Naval  Research.  GAV  is  a  recipient  of  a  National  Science  Foundation  Presidential  Young 
Investigator  Award,  a  David  and  Lucile  Packard  Fellowship  in  Science  and  Engineering,  an 
Alfred  P.  Sloan  Foundation  Research  Fellowship,  and  a  Dreyfus  Foundation  New  Faculty 
Award. 


-28- 


Appendix:  The  Renormalized  Centroid-Constrained  Propagator  in  Multidi¬ 
mensional  Space 

It  is  first  assumed  that  the  reference  centroid-constrained  propagator  [cf.  Eq.  (3.2)] 
for  a  multidimensional  system  is  a  diagonal  matrix,  that  is.  A{t)  =  where  the 

uppercase  letters  here  stand  for  matrices  or  vectors.  Because  there  are  mixed-index  partial 
derivatives  at  the  vertices,  the  multidimensional  centroid-constrained  propagator  Cdr.qc) 
matrix  [cf.  Eq.  (2.1)]  is  not  necessarily  diagonal.  (Note  here  that  Qc  is  now  the  multidimen¬ 
sional  cnutroid  variable.)  The  first  coirection  to  consider  ia  the  chain  summation  in  Eq. 

(3.36)  which  gives  rise  to  the  optimized  centroid-constreiined  correlation  function  matrix 

Ccir,qc)  =  A-i3{{AD):{DCc)AV)c,  ,  (Al) 

where  D  is  the  partial  derivative  vector  and  AP  is  the  scalar  multidimensional  po¬ 

tential  difference.  The  notation  (•  •  •)c^  here  denotes  a  multidimensional  Gaussian  average 
with  width  factors  taken  from  the  matrix  of  optimized  Cc(0,  qd  values.  Upon  substituting 
the  explicit  form  of  A  in  Eq.  (A.l),  one  obtains  the  multidimensional  matrbc  analog  of  Eq. 

(3.37)  as 


which  leads  to  the  definition  of  the  optinuzed  frequency  tensor, 

=  {dtxd^y)cc  ■  (^-3) 

This  is  the  same  self-consistent  equation  obtained  previously  [see,  e.g.,  Eq.  (2.48)  of  Ref. 
10].  The  two-line-loop  correction  in  Eq.  (3.38)  is  more  complicated,  but  as  long  as  one  uses 
the  tensor  prescription  it  is  always  possible  to  generalize  the  one-dimensional  equations  to 
multidimensional  space. 
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Table  I:  The  average  value  {q~)  for  the  potential  in  Eq.  (4.1)  evaluated  with  the  centroid- 
based  formalism  in  Eq.  (2.18)  and  with  various  levels  of  approximation  for  the  centroid 
density.^'  The  exact  results  were  obtained  by  path  integral  Monte  Carlo. 


d 

{Q^)  exact 

{q^)i 

{q^)2 

61 

60 

5 

1.378 

1.361 

1.371 

1.3 

0.5 

7 

1.382 

1.352 

1.372 

2.2 

0.7 

9 

1.384 

1.344 

1.372 

2.9 

0.9 

“'The  results  {q~)i  based  on  the  optimized  LHO  reference  potential  approximation  for 
the  centroid  density  in  Eqs.  (3.22)-(3.26),  while  {q^)2  are  the  results  including  the  higher- 
order  corrections,  Eqs.  (3.28)  and  (3.29).  The  quantities  (5i  and  62  are  the  percentage 
errors  of  {q~)i  and  (g^)2,  respectively,  compared  with  the  numerically  exact  result. 


-32- 


Table  II:  Quantum  Correction  Factors“^  from  Eq.  (4.2)  for  the  Eckart  barrier.^' 


u 

El 

r2 

Ea/c 

6 

4.4 

4.4 

4.4 

8 

15.0 

17.0 

17.0 

10 

73.0 

110.6 

105.0 

12 

514.0 

1278.0 

1240.0 

“’The  quantum  corrections  E i  are  based  on  the  optimized  LHO  reference  potential  approx¬ 
imation  for  the  centroid  density  in  Eqs.  (3.22)'(3.26),  while  r2  are  the  results  including 
the  higher-order  terms  in  Eqs.  (3.31).  The  quantum  correction  Ea/c  is  the  path  integral 
Monte  Carlo  result  reported  in  Ref.  12a. 

^’The  Eckart  barrier  potential  is  given  by  V (q)  =  Vosech"^{q/aQ),  with  the  parameter  values 
'IirVo/huf,  =  12.0  and  u  =  dhut,  in  the  present  calculations,  and  Ub  is  the  magnitude  of  the 
classical  barrier  frequency. 
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Figure  Captions 


Fig.  1:  A  schematic  diagram  depicting  a  discretized  Feynman  path  q{T)  within  the  imaginary 
time  interval  0  <  r  <  Tip.  The  isomorphic  classical  quasiparticles  are  shown  by  the 
dark  circles  which  form  a  “ring  polymer” .  Each  quasiparticle  “interacts'’  with  its  two 
nearest- neighbors  through  effective  harmonic  forces  and  feels  the  external  potential 
through  the  term  V{qi)/P  [cf.  Eq.  (1-2)].  The  centroid  variable  defined  in  Eqs. 
(1.5)  and  (1.6)  is  also  shown. 

Fig.  2:  A  schematic  diagram  depicting  the  calculation  of  an  operator  average  (A)  in  the 
discrtetized  Feynman  path  integral  picture.  The  operator  A[q{T)]  can  be  evaluated 
anywhere  on  the  path  integral  ring  0  <  t  <  h3  due  to  the  cyclic  invariance  of  the 
trace  operation.  The  centroid  variable,  which  is  off  the  ring,  is  also  shown. 

Fig.  3:  A  schematic  diagram  depicting  the  calculation  of  the  imaginary  time  correlation  func¬ 
tion  (A(r)S(O))  in  the  discretized  Feynman  path  integral  picture.  The  operators 
A[q{r')]  and  B[q{T")]  can  be  evaluated  at  any  two  points  on  the  path  integral  ring 
subject  to  the  constraint  0  <  r  =  r'  —  r"  <  hp.  The  centroid  variable,  which  is  off 
the  ring,  is  also  shown. 

A  plot  of  the  imaginary  time  correlation  function  (9^(T)q^(0))  for  the  non-quadratic 
potential  described  in  Sec.  IV  [Eq.  (4.1)].  The  correlation  function  is  plotted  as  a 
function  of  the  dimensionless  variable  u  =  r/hP  with  P  =  b.  The  solid  circles  show 
the  numerically  exact  results  while  the  solid  line  is  for  the  optimized  LHO  theory  in 
Eqs.  (3.22)  -  (3.26)  with  the  centroid-based  formulation  of  the  correlation  fimction  in 
Eq  (2.18). 

A  plot  of  the  imaginary  time  correlation  function  (q^(r)g^(O))  for  the  non-quadratic 
potential  described  in  Sec.  IV  [Eq.  (4.1)].  The  correlation  function  ie  plotted  as  a 
function  of  the  dimensionless  variable  u  =■  r/HP  with  /?  =  5.  The  solid  circles  show 
the  niimerically  exact  results  while  the  solid  line  is  for  the  optimized  LHO  theory  in 


Fig.  4: 


Fig.  5: 
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Eqs.  (3.22)  -  (3.26)  with  the  centroid-bcised  formulation  of  the  correlation  function  in 
Eq  (2.18). 

Fig.  6:  A  plot  of  the  centroid-constrained  correlation  function  Cc(t,  Qc)  defined  in  Eq.  (2.1) 
for  the  non-quadratic  potential  described  in  Sec.  IV  [Eq.  (4.1)]  The  numerically  ex¬ 
act  results  are  shown  by  the  solid  circles.The  correlation  function  obtained  from  the 
optimized  LHO  approximation  in  Eq.  (3.37)  is  shown  by  the  dashed  line,  while  the 
solid  line  shows  the  results  obtained  by  including  the  two-line-loop  correction  from 
Eq.  (3.40).  The  correlation  frmctions  are  plotted  as  a  function  of  u  =  r/h(3  and  for 
/3  =  10  and  Qc  =  0.0. 

Fig.  7;  A  plot  of  the  imaginary  time  position  correlation  function  C(t)  in  Eq.  (2.3)  for  the 
non-quadratic  potential  described  in  Sec.  IV  [Eq.  (4.1)].  The  numerically  exact  results 
are  shown  by  the  sohd  circles.  The  correlation  fxmction  obtained  using  the  centroid- 
constrained  optimized  LHO  approximation  in  Eq.  (3.37)  is  shown  by  the  dashed  line, 
while  the  solid  Une  shows  the  results  obtained  by  including  the  two- line-loop  correction 
from  Eq.  (3.40).  In  the  two  latter  cases,  the  correlation  function  C(r)  was  obtained 
by  averaging  the  appropriate  analytical  centroid-constrained  correlation  function  over 
the  numerically  determined  centroid  density.  The  results  are  plotted  as  a  function  of 
u  =  r/hP  and  for  0  =  5. 
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Centroid-Constrained  Imaginary  Time  Correlation  Function 
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